1 Notebook settings and libraries

set.seed(123)
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
knitr::opts_chunk$set(fig.width = 10, fig.height = 10, fig.dpi = 600)
require(EnhancedVolcano)
require(DESeq2)
require(limma)
require(IHW)
require(tidyverse)
require(dplyr)
require(ggplot2)
require(Hmisc)
require(org.Hs.eg.db)
require(clusterProfiler)

2 Global Variables

SIG_THR <- 0.05
LFC_THR <- 1
MIN_COUNT <- 10
MIN_SUMMED_ROW_COUNT <- 3

3 Functions

# Source: https://rdrr.io/github/ttecon/ttr/src/R/geom_split_violin.R
#
#' Split Violin plot
#'
#' This function is to plot the split violin graph.
#' It is built on [ggplot2::geom_violin()],
#' and adopted from jan-glx's stackoverflow answer \url{https://stackoverflow.com/a/45614547}
#' A violin plot is a compact display of a continuous distribution.
#' It is a blend of [ggplot2::geom_boxplot()] and [ggplot2::geom_density()]:
#' A violin plot is a mirrored density plot displayed in the same way as a boxplot.
#' A split violin plot allows for odd violins to be plotted to the left,
#' and even violins to be plotted to the right
#'
#' @param mapping Set of aesthetic mappings created by `aes()` or `aes_()`.
#' If specified and `inherit.aes = TRUE` (the default),
#' it is combined with the default mapping at the top level of the plot.
#' You must supply mapping if there is no plot mapping.
#' @param data The data to be displayed in this layer. There are three options:
#'
#' If `NULL`, the default, the data is inherited from the plot data
#' as specified in the call to [ggplot2::ggplot()].
#'
#' A `data.frame`, or other object, will override the plot data.
#' All objects will be fortified to produce a data frame.
#' See [ggplot2::fortify()] for which variables will be created.
#'
#' A function will be called with a single argument, the plot data.
#' The return value must be a `data.frame`, and will be used as the layer data.
#' A function can be created from a formula (e.g. `~ head(.x, 10)`).
#' @param stat Use to override the default connection between `geom_density` and `stat_density`.
#' @param position Position adjustment, either as a string,
#' or the result of a call to a position adjustment function.
#' @param ... Other arguments passed on to `layer()`.
#' These are often aesthetics, used to set an aesthetic to a fixed value,
#'  like `colour = "red"` or `size = 3`.
#' They may also be parameters to the paired geom/stat.
#' @param draw_quantiles If `not(NULL)` (default), draw horizontal lines
#' at the given quantiles of the density estimate.
#' @param trim If `TRUE `(default), trim the tails of the violins to the range of the data.
#' If `FALSE`, don't trim the tails.
#' @param scale If `"area"` (default), all violins have the same area (before trimming the tails).
#' If `"count"`, areas are scaled proportionally to the number of observations.
#' If `"width"`, all violins have the same maximum width.
#'
#' @param na.rm If `FALSE`, the default, missing values are removed with a warning.
#' If `TRUE`, missing values are silently removed.
#' @param show.legend logical. Should this layer be included in the legends?
#' `NA`, the default, includes if any aesthetics are mapped.
#' `FALSE` never includes, and `TRUE` always includes.
#' It can also be a named logical vector to finely select the aesthetics to display.
#' @param inherit.aes If `FALSE`, overrides the default aesthetics,
#' rather than combining with them.
#' This is most useful for helper functions that define both data and aesthetics
#' and shouldn't inherit behaviour from the default plot specification, e.g. [ggplot2::borders()].
#'
#' @export
geom_split_violin <- function(mapping = NULL,
                              data = NULL,
                              stat = "ydensity",
                              position = "identity",
                              ...,
                              draw_quantiles = NULL,
                              trim = TRUE,
                              scale = "area",
                              na.rm = FALSE,
                              show.legend = NA,
                              inherit.aes = TRUE) {
  ggplot2::layer(
    data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin,
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(
      trim = trim, scale = scale, draw_quantiles = draw_quantiles,
      na.rm = na.rm, ...
    )
  )
}

#' @format NULL
#' @usage NULL
GeomSplitViolin <- ggplot2:::ggproto("GeomSplitViolin",
  ggplot2::GeomViolin,
  draw_group = function(self, data, ..., draw_quantiles = NULL) {
    data <- transform(data, xminv = x - violinwidth * (x - xmin), 
                      xmaxv = x + violinwidth * (xmax - x))
    grp <- data[1, "group"]
    newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv), 
                             if (grp %% 2 == 1) y else -y)
    newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
    newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
    if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
      stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
        1))
      quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
      aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), 
                         drop = FALSE]
      aesthetics$alpha <- rep(1, nrow(quantiles))
      both <- cbind(quantiles, aesthetics)
      quantile_grob <- GeomPath$draw_panel(both, ...)
      ggplot2:::ggname("geom_split_violin", 
                       grid::grobTree(GeomPolygon$draw_panel(newdata, ...), 
                                      quantile_grob))
    }
    else {
      ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
    }
  }
)

4 Differential gene expression analysis in HGSOC samples

4.1 Data preparation

Feature counts were obtained using Rsubrea::featureCounts(). Strandness was inferred using salmon utilized in istrand script. All samples were inward inversely stranded libraries.

# load count matrix
counts <- read.csv("workdir/feature_counts/feature_counts.csv", 
                   row.names = 1)
cols <- colnames(counts)
cols <- gsub("_sorted.bam", "", cols)
cols <- gsub("X", "", cols)
colnames(counts) <- cols

# change entrezIDs to gene symbols if available
id.to.symb <- bitr(rownames(counts), 
                   fromType = "ENTREZID",
                   toType = c("SYMBOL"),
                   OrgDb = org.Hs.eg.db,
                   drop = FALSE)

gene.symbols <- ifelse(is.na(id.to.symb$SYMBOL), 
                       id.to.symb$ENTREZID, 
                       id.to.symb$SYMBOL)

rownames(counts) <- gene.symbols
counts[1:5, 1:5]
df <- read.csv("workdir/clinical.csv")

df$OS_median <- ifelse(df$OS > median(df$OS), "above", "below") %>% 
  as.factor() %>% relevel(ref = "above")

df$DFS_median <- ifelse(df$DFS > median(df$DFS), "above", "below") %>% 
  as.factor() %>% relevel(ref = "above")

df$CR <- df$CR %>% as.factor

str(df)
## 'data.frame':    165 obs. of  17 variables:
##  $ ID           : chr  "8" "11" "30" "37" ...
##  $ therapy      : chr  "PC" "PC" "PC" "PC" ...
##  $ age          : int  55 65 34 48 46 66 59 64 45 53 ...
##  $ CR           : Factor w/ 2 levels "no","yes": 1 2 2 2 2 1 1 1 2 2 ...
##  $ PS           : chr  "no" "yes" "yes" "yes" ...
##  $ OS           : int  104 1193 467 2320 3750 385 790 1138 1847 2351 ...
##  $ DFS          : int  0 393 226 815 1468 0 0 0 1429 528 ...
##  $ relapse      : chr  NA "yes" "yes" "yes" ...
##  $ death        : chr  "yes" "yes" "yes" "yes" ...
##  $ TP53_mm      : chr  "yes" "yes" "no" "no" ...
##  $ p53_acc      : chr  "yes" "yes" "no" "no" ...
##  $ type         : chr  "HGSOC" "HGSOC" "LGSOC" "LGSOC" ...
##  $ microinvasion: chr  NA NA NA NA ...
##  $ implant      : chr  NA NA NA NA ...
##  $ primary_tumor: chr  NA NA NA NA ...
##  $ OS_median    : Factor w/ 2 levels "above","below": 2 2 2 1 1 2 2 2 1 1 ...
##  $ DFS_median   : Factor w/ 2 levels "above","below": 2 1 2 1 1 2 2 2 1 1 ...

4.2 OS median below vs above

Because of the difference in expression profiles between HGSOC and other types and the fact that HGSOC is the largest sample group this group was selected for subsequent analysis. Number of HGSOC samples was 140.

dff <- df %>% filter(type == "HGSOC")
nrow(dff)
## [1] 140

Firstly, patients with OS above median were compared to those with OS below median. Low count rows (genes) were filtered out to avoid noise. Exact experiment design is shown below.

dds <- DESeqDataSetFromMatrix(countData = counts[dff$ID],
                              colData = dff,
                              design = ~OS_median)

keep <- rowSums(counts(dds) >= MIN_COUNT) >= MIN_SUMMED_ROW_COUNT
dds <- dds[keep, ]

dge <- DESeq(dds)
res <- results(dge, filterFun = ihw)
head(res, 0)
## log2 fold change (MLE): OS median below vs above 
## Wald test p-value: OS median below vs above 
## DataFrame with 0 rows and 7 columns

Number of significantly differentiating genes and volcano plot are show below.

signif <- res %>% as.data.frame %>% filter(abs(log2FoldChange) > LFC_THR,
                                           padj < SIG_THR)
write.table(signif, file = "workdir/DGEA/prog.tsv", quote = FALSE)
length(signif %>% row.names)
## [1] 219
EnhancedVolcano(res, 
                lab = rownames(res), 
                x = 'log2FoldChange', 
                y = 'padj',
                xlab = "LFC",
                ylab = "FDR",
                pCutoff = SIG_THR,
                FCcutoff = LFC_THR,
                legendLabels = c("NS", "LFC", "FDR", "FDR and LFC"))

resLFC <- lfcShrink(dge, coef = 2, type = "apeglm")
DESeq2::plotMA(resLFC, ylim = c(-5, 5))

Top 10 significant genes were selected and the differences between 2 groups are plotted below according to (a) log count (b) log SF normalized count (c) VST normalized count.

top_df <- signif %>% arrange(padj) %>% head(10) 
top_df
dds <- estimateSizeFactors(dds)
norm.counts <- counts(dds, normalized = TRUE)
vst.counts <- dds %>% vst %>% assay

top_genes <- top_df %>% rownames

top_genes_data <- data.frame()

for (gene in top_genes) {
  d <- plotCounts(dds, gene=gene, intgroup="OS_median", returnData=TRUE)
  d$ID <- rownames(d)
  d$gene <- gene
  new_df <- data.frame(gene=d$gene, ID=d$ID, 
                       OS_median=d$OS_median, 
                       count=d$count)
  top_genes_data <- rbind(top_genes_data, new_df)
}

top_genes_data$norm_count <- (norm.counts[unique(top_genes_data$gene), ] 
                              %>% t)[dff$ID,] %>% as.vector
top_genes_data$vst_count <- (vst.counts[unique(top_genes_data$gene), ] 
                             %>% t)[dff$ID,] %>% as.vector
ggplot(top_genes_data, aes(x=gene, y=log(count), fill = OS_median)) +
  geom_split_violin(trim = 0) +
  geom_boxplot(alpha = 0.3, width = 0.2, outliers = FALSE) +
  theme_classic() +
  ylab("log count") +
  guides(fill = guide_legend(title = "OS median")) +
  theme(axis.line.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.x = element_blank())

ggplot(top_genes_data, aes(x=gene, y=log(norm_count), fill = OS_median)) +
  geom_split_violin(trim = 0) +
  geom_boxplot(alpha = 0.3, width = 0.2, outliers = FALSE) +
  theme_classic() +
  ylab("log SF normalized count") +
  guides(fill = guide_legend(title = "OS median")) +
  theme(axis.line.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.x = element_blank())

ggplot(top_genes_data, aes(x=gene, y=vst_count, fill = OS_median)) +
  geom_split_violin(trim = 0) +
  geom_boxplot(alpha = 0.3, width = 0.2, outliers = FALSE) +
  theme_classic() +
  guides(fill = guide_legend(title = "OS median")) +
  ylab("VST normalized count") +
  theme(axis.line.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.x = element_blank())